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Abstract 

O . 

The Navier-Stokes equations in a finite cylinder are written in terms of poloidal 
^ \ and toroidal potentials in order to impose incompressibility. Regularity of the so- 
^ \ lutions is ensured in several ways: First, the potentials are represented using a 
(-H \ spectral basis which is analytic at the cylindrical axis. Second, the non-physical 
' discontinuous boundary conditions at the cylindrical corners are smoothed using 

a polynomial approximation to a steep exponential profile. Third, the nonlinear 
term is evaluated in such a way as to eliminate singularities. The resulting pseudo- 
spectral code is tested using exact polynomial solutions and the spectral conver- 
. ^ ! eence of the coefficients is demonstrated. Our solutions are shown to agree with 
00 , exact polynomial solutions and with previous axisymmetric calculations of vortex 
^ \ breakdown and of nonaxisymmetric calculations of onset of helical spirals. Paral- 
^ ' lelization by azimuthal wavenumber is shown to be highly effective. 

o 

^ \ Key words: cylindrical coordinates, coordinate singularities, pseudo-spectral, 
recursion relations, radial basis, vortex breakdown, spectral convergence. 
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1 Introduction 



The von Karman flow owes its name to T. von Karman [1] who in 1921 flrst studied 
the flow in the semi-inflnite domain bounded by a single rotating disk using a simi- 
larity transformation. In 1951, Batchelor [2] extended the problem to the flow confined 
between two infinite rotating disks. For rotating disks of finite radius, the configura- 
tion can be described by three control parameters: the ratio s of the angular velocity 
of the two disks, the height-to-radius ratio h and a Reynolds number Re based on the 
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radius and the azimuthal velocity of one of the disks. The variation of these three pa- 
rameters {Re,s,h) has proved to yield a rich variety of qualitatively different accessible 
flows, even before the onset of turbulence. The symmetries influence the transitions 
that the flow can undergo. This configuration is extensively studied in the context of 
transition to complex and turbulent flows. All of these properties explain why the von 
Karman flow is increasingly considered as one of the classical hydrodynamic configu- 
rations and why the scientific community is interested in further exploring its complex 
behavior. 

The first numerical studies were necessarily devoted to axisymmetric flows and their 
stability. In the rotor-stator configuration (s — 0), vortex breakdown forming character- 
istic recirculation bubbles was observed by number of authors (Lugt and Abboud [3], 
Daube and Sorensen [4], Lopez [5], Daube [6]). This now well-documented configura- 
tion became a benchmark for testing axisymmetric codes. Following Lopez and Shen 
[7], Speetjens [8] and a number of other authors, we will validate our method in the 
axisymmetric configuration by reproducing the stationary state at Re = 1800 and, for 
Re = 2800, the non-stationary, oscillating flow, for which we will compare the bifurca- 
tion threshold and the oscillation frequency against previous findings. 

Increasing computational power has made it possible to study three-dimensional in- 
stabilities. Breaking of axisymmetry has been the subject of several studies, of which 
we mention these of Gauthier et al. [9], Gelfgat et al. [10], Blackburn and Lopez [11], 
Lopez et al. [12] and Nore et al. [13]. Three-dimensional instability precedes axisym- 
metric instability for h <1.6 and h > 2.8. As the test problem for validating our code in 
three dimensions, we have selected a configuration with (s = — = 3.5, Re = 2150), 
where the perturbation of the axisymmetric flow takes the form of a helical spiral. In 
section 3.5 we compare our results with those of Lopez et al. [14] and Gelfgat et al. [10]. 

Interesting phenomena can also be observed in the turbulent regime. Turbulence and 
large-scale structures may coexist. In experiments in a turbulent counter-rotating con- 
figuration, Marie [15] and Ravelet et al. [16] discovered that a two-cell mean flow with 
a shear layer at the cylinder mid-plane undergoes switching to a one-cell mean flow 
whose shear layer is adjacent to the less rapidly rotating disk. This transition can be 
observed at Reynolds numbers which are numerically accessible. 

A comprehensive classification of the solutions for different values of the parame- 
ters (s,h,Re) is beyond the scope of this work. Our main purpose here is to develop 
a mathematical and algorithmic tool which can be applied to von Karman flow and 
to rotating turbulence, and which can be extended to the magnetohydrodynamic con- 
figuration of the VKS experiment [17]. The major component of our algorithm is the 
poloidal-toroidal decomposition [18, 19], which insures incompressibility by construc- 
tion, at the price of increasing the order of the governing equations. When applied to 
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the Navier-Stokes equation in a finite cylinder, the resulting system has boundary con- 
ditions which are coupled and of high order. In a companion article [20], we showed 
that this system could be reduced to the solution of a set of nested Helmholtz and Pois- 
son problems with uncoupled Dirichlet boundary conditions, whose solutions could 
be superposed via the influence matrix technique. The purpose of the present article 
is to describe our method for solving these elliptic problems using a spectral represen- 
tation which exploits azimuthal symmetry of the system and which is regular at the 
cylindrical axis, and to demonstrate the validity of the resulting hydrodynamic code. 

More specifically, it was shown by Marques and co-workers [18, 19] that the Navier- 
Stokes equations 

dtu + (u • V)u = Re"^Au - Vp (1.1a) 
V • u = (1.1b) 

in a finite cylinder with boundary conditions 

u = rco±e0 at z = ±^, (l-2a) 
u = at r = l, (1.2b) 

and with toroidal and poloidal potentials defined by 

u = V X (ipez) + V X V X {(pez) (1.3) 

are equivalent to the two scalar equations 

(dt - Re-'^A)Ah^ = S^p= ez • V X (u • V)u (1.4a) 
{dt - Re~'^A)AAh(p = S^ = -ez • V x V x (u • V)u (1.4b) 

where = jdyrdr + ^3^, with boundary conditions 
1 1 

-de^ + ^rz(p = = A}i<p = (p = drzA^xp — -dgAA^cf) = at r = 1 (1.5a) 

^^0 at r = (1.5b) 

1 ft 

Afif + -dr (^co±r^^ — ^z^h<P = ^h<P = at z = ±- (l-5c) 

Our article [20] was devoted to showing how the problem (1.4)-(1.5) can in turn be re- 
duced to a sequence of five nested parabolic and elliptic problems, each with Dirichlet 
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boundary conditions: 



dt - ftp = Srp f^\^^i = (rf fxp\^^^h = --drico±r^) (1.6a) 

\^=U drXp\r=l=0 tp\r=0 = (1.6b) 

(dt-Re-^A^g = S^ g\r=i = (rg g\^^^=(r^ (1.6c) 

AU^g Mr=i = M,=±|=0 (1.6d) 

Ah(p^U (p\r=i^O (l-6e) 



for the two potentials ip, (p and three intermediate fields g, f^p and The influence 
matrix technique [21], a generalization of the usual separation into particular and ho- 
mogeneous solutions, is used to determine boundary values df, dg, such that the 
boundary conditions present in (1.5) but not in (1.6) are satisfied, i.e. such that 

1 (1.7a) 
1 (1.7b) 
±^ (l-7c) 

In this article, we describe the numerical implementation of this algorithm. We first 
present the spatial discretization of the fields, using a set of basis functions [22] which is 
regular at the cylindrical axis r — 0, and regularizing the discontinuous boundary con- 
ditions at the corners r — l,z — ±h/2. We then explain the methods we have used for 
solving equations (1.6), in particular for stably and economically solving the Helmholtz 
problems resulting from time discretization and for evaluating the nonlinear terms. Fi- 
nally, we describe the validation of the implementation, comparing results from our 
code to an analytic polynomial solution and to previously published two- and three- 
dimensional test cases. 



^rzftp - -^eg =0 at r = 
1 

-dexp + drz(p =0 at r = 

3z/^ =0 at z = 



2 Spatial and temporal discretization 

The spectral discretization that we use is 

I M I I M I 

/(r,0,z)= E r(r,z).-^= L L L me'^'Qnirml-) (2.1) 

m=-|f I m=-\^ \ k=0 n=\m\ V « / 

In (2.1), we do not introduce new notation for Fourier coefficients, or for coefficients 
in the 3D tensor-product basis, using instead the number and type of superscripts and 
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subscripts to distinguish between functions in physical space and spectral space coef- 
jficients. 

The basis functions in the azimuthal and axial directions are standard [23, 24]: Fourier 
modes e''"^ and Chebyshev polynomials T]^{2z/h), respectively. In cylindrical geome- 
tries, it is the radial direction which is most problematic and on which we will focus. To 
represent this direction, we use the basis functions (r) developed by Matsushima 
and Marcus [22]. In sections 2.1 and 2.2, we will discuss the means by which we impose 
regularity at the origin r = and at the comers r = 1. 

2.1 Regular basis of radial polynomials 

A function / on the disk is analytic at the origin if the radial dependence of the Fourier 
coefficient /'"(r) multiplying the Fourier mode e^^^ is of the form 

oo 

f^{r,z)= C{z)r^ = ccZ{z)r^ + ccl^2{^)r"'+^ + ---^r^p{r^) (2.2) 

n=\m\ 
n+m even 

where pis a polynomial. Examples of functions which violate (2.2) are given in figures 

1. 2 and 3 and compared with regular functions obeying (2.2) on the right of each figure. 

Various approaches used in spectral methods to represent functions in polar coordi- 
nates are surveyed by Boyd [25, 26] and by Canuto et al. [24]. A common practice 
[14, 21, 27, 28] has been to impose some degree of continuity, such as C^, but not com- 
plete analyticity C°°. Although basis functions which are not analytic at the origin gen- 
erally do not pollute the fields, retaining such functions wastes CPU time and memory 
which could be better used to increase resolution. 

The condition (2.2) is stated in terms of monomials, the use of which is excluded be- 
cause of their poor numerical properties. The polynomial basis developed by Mat- 
sushima and Marcus [22] respects these conditions and yet is numerically well-conditioned. 
These polynomials QJf (a,jS;r) are solutions to the singular Sturm-Liouville equation: 

f (1 - r^V'"^ d / / ?\'' P,d\ \m\(\m\ + B-l) , ^ ^ ^m. ^ n ^ 

y J Tr [[^ - n '^TrJ ~ r2 ' + n{n + lcc + ^-l)^QZ{a,^-r)=Q, 

(2.3) 

defined over r G [0, 1] . In (2.3), < | m | < n, a; G [0, 1] and /3 is a positive integer. With the 
special choice a = jS = 1, Q^(l,l;r) are related to Legendre and shifted Jacobi polyno- 
mials used by Leonard and Wray [29]; similar functions were also derived by Verkley 
[30]. The functions QJf (a,jS;r) are complete and orthogonal over [0,1] with respect to 
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Fig. 1. Coordinate singularity effects: parity mismatch. Left: f{r,6) = r. Right: f{r,9) = r^. 
f{r,e)=r^ + ^cos{4e) f{r,e)^r^ + ^cos{49) f{r, 6) = + ^cos{Ae) 




Fig. 2. Coordinate singularity effects. Left: discontinuity of value. Middle: discontinuity of 
Laplacian. Right: regular function. 
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Fig. 3. Clustering effect - contours of f'"{r,9) = 0.5. Left: f"'{r,e) =r^ + 32r^cos{me), which is 
not smooth at for m > 2. Right: f" (r, 9) = r^ + 32r'" cos{m6), which is smooth at 0. 
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the inner product: 

I' Q;r(«.iS;r)Q-,(«,)S;r)^^-^^ = (2.4) 
The order polynomials Q^{ix,^;r) have the following explicit expression: 

Jn ^^,P,r) _ + - p + 1) r + v) r (^^) 

(2.5) 

but they, as well as the normalizing coefficients 7^ ( a, /3) , can be calculated in O (n — \m\) 
operations using recursion relations given by Matsushima and Marcus [22]. Recursion 
relations also exist for the operators 

\r^r,r^,ir^rf - m^,{r^rf + Ar^j (2.6) 

expressed in the QJf polynomial basis, meaning that for any of the operators H in (2.6), 
there exist banded matrices L and R such that H = R~^L. Thus, 

Hf^g<^Lf^Rg (2.7) 

reducing the time for multiplication by H or from quadratic to linear in the num- 
ber of radial modes or gridpoints. The existence of recursion relations is a general prop- 
erty of differential operators represented in polynomial bases, for reasons explained by 
Tuckerman [31]. Recursion relations will be further discussed in section 2.4. 

The radial function f'^ (r) associated with Fourier mode m and its coefficients in the 
polynomial basis are related by the transform pair: 

00 N 

rW= L fnOZir)^ E /-Q-(r) (2.8a) 



n=\m\ n=\m\ 
n+meven n+meven 



fl N-1 

= / rfrw;(r)r (r)Q-(r)/;- « ^ Wiriri)Q^iri)/C (2.8b) 



In (2.8), the order of the polynomial expression is N = 2N — 2 or N = 2N — 1 accord- 
ing to whether m is even or odd. The collocation points {r/} for Gaussian quadrature 
are computed numerically as the roots of the first neglected m = polynomial Q^^2 
if the boundary points are to be excluded; otherwise, they are the roots of a slightly 
more complicated expression [22]. Once the {rj} are determined, the weights {w/} are 
computed by recursion relations [22]. 

Equations (2.8) specify Nm = N — [m/2] coefficients from values at N quadrature points 
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via a rectangular matrix. Since the basis is orthonormal, the inverse transformation is 
obtained from the transpose of this rectangular matrix. The smaller size of the spectral 
representation compared to the grid representation is a consequence of the fact that the 
functions in (2.8) are not arbitrary functions of r, but belong to the restricted subspace 
of functions obeying the regularity conditions (2.2). 

2.2 Regularization of the comers 

The boundary conditions on uq stated in (1.2) are 



The equations have been nondimensionalized using the radius as the unit of length 
and the inverse angular velocity 1/co- (with a;_ > and |a; + | < a;-) as the unit of 
time. Therefore, cv- = 1 and — 1 < s = co+ < 1, in particular cv+ is 0, +1 or —1, for the 
rotor-stator, exactly corotating, or exactly counter-rotating configurations, respectively. 
(It is also possible to set the velocity on the cylinder to some other constant value 
instead of 0; for example to simulate the flow in a rotating cylinder, we would set 
Me|r=i(0,z) = C0+ = co^ = 1. Here, we focus instead on the difficulties in implementing 
the discontinuous boundary conditions (2.9).) 

Boundary conditions (2.9) are discontinuous at the comer points r — 1, z — —h/2 and 
possibly r — 1, z — +h/2: one or both disks rotate while the lateral boundary remains 
fixed. Mathematically, a PDE with a finite number of singular points can have a solu- 
tion which is smooth except at these points. However, spectral methods then do not 
converge exponentially because series of smooth functions cannot converge uniformly 
to a discontinuous solution. If nothing is done to prevent it, the Gibbs phenomenon 
will lead to spurious oscillations which propagate into the whole domain from the 
neighborhood of the singularity. For finite difference methods, the discontinuity will 
affect only the neighborhood of the singular point, on the order of the grid interval, and 
therefore does not pose a severe problem. Finite volume methods have a local integral 
formulation and so the discontinuity presents an even less serious problem. The filter- 
ing intrinsic to local methods is, however, intrinsically related to the high numerical 
diffusion which in turn makes local methods less precise. 

In some cases, even local methods do not sufficiently filter singularities. Georgiou et al. 
[32] discuss the issue of spurious oscillations in the context of finite element methods. 
In the solutocapillary problem studied by Martin-Witkowski and Walker [33], the au- 
thors were required to explicitly filter the solution to achieve acceptable convergence 



ue{e,z) =0 
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(2.9a) 
(2.9b) 
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even in a finite difference calculation. 



Spectral methods must always explicitly filter strong singularities like (2.9a)-(2.9b). We 
have chosen to do this by approximating the discontinuous function at the boundary 
by a steep but smooth profile. This procedure can be justified by arguing that we are not 
interested in finding the solution to the singular problem. In any real experiment, the 
boundary conditions are not discontinuous: a small gap must necessarily exist between 
the rotating disks and the stationary lateral boundary. In our algorithm, we replace the 
constant angular velocities in (2.9a) by continuous functions, as illustrated by figure 
4a. Two possibilities are: 

U0{r,6) —r (^1 — e"^^ cci± at z 
U0{r,6) =r{l -r^)co± at z 

where }i is an arbitrary but large even integer (e.g. ji = 10), and S, iv± are constants. 
In order to be represented in a radial polynomial basis, the exponential regularization 
(2.10a) must be approximated by a polynomial, so that both expressions above are 
effectively polynomials. The steepness of the profiles are adjusted by varying 3 or ji, 
whose possible values are limited by the radial polynomial order N. 

The advantage of exponential regularization is that, for a given N, a steeper profile 
can be achieved by the polynomial approximation to (2.10a) than by the polynomial 
(2.10b), as can be seen in figure 5a. In this way, the deviation from the idealized profile 
is minimized, without unduly increasing the radial resolution N over that required to 
resolve the field in the interior. The polynomial approximation to (2.10a) differs from 1 
by more than 10% only over the small range 1 — 2^ < r < 1. 

Regularization of the boundary condition imposes a lower bound on the spectral res- 
olution - the spectral approximation must be able to represent the regularization pro- 
files smoothly. In fact, Lopez and Shen [7] observed that the actual resolution should 
be approximately twice the minimal resolution sufficient for representing the regular- 
ization profiles, due to generation of higher wavenumbers by the nonlinear terms. As 
was shown by Lopez and Shen [7], for comparing with results obtained using different 
methods and for benchmark purposes it is sufficient to use 3 ^ 0.005. In practice, we 
use 0.005 <S < 0.05, illustrated in figure 5b. 

Care must be taken in approximating (2.10a) by a polynomial expression, in order to 
satisfy all of the conditions that we require of U0 and of tp. The procedure we use is 
as follows. We evaluate (2.10a) on the collocation points. Since ue is an odd function 
of r (see (2.30)), we apply the transform (2.8) using the (odd) pol5momials associated 
with the m — 1 Fourier mode. The basis of odd radial functions insures that this ap- 
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(2.10a) 
(2.10b) 
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proximation to on the bounding disks is zero at r = 0, while the use of Gauss-Radau 
collocation points which include the cylinder boundary insures that that it is also zero 
air — 1. This approximation to ug is integrated over r to obtain an approximation to tp, 
with the integration constant chosen in order to satisfy — air — 0. 




Fig. 4. Regularized profiles used for elimination of the discontinuous boundary conditions at 
the cylinder corners, a) regularization on upper and bottom disks using (2.10). b) regularization 
on lateral bounding cylinder using (2.11). 




0.2 0.4 0.6 0.8 1 ■^^ 0.2 0.4 0.6 0.8 1 



Fig. 5. a) Comparison between regular profiles represented by polynomials of order r^. Solid 
curve: uqI^^^h (r) = r (l — r^). Dotted curve: uqI^^^h (r) ^r(l — Inset illustrates result 

b) Comparison between exponential profiles ue\^^^h (r) =r(l— for different S. 
Another choice is to apply a filter to the lateral boundary, replacing (2.9b) by 

ueiz)=cv+e-(^~T)/^ + cv-e-(^+^)/^ at r = l (2.11) 

while keeping the uniform angular velocity profiles (2.9a) on the disks unchanged; see 
figure 4b. This kind of regularization is similar to that of Lopez and Shen [7]. 

10 



Other, quite different, approaches to the treatment of singularities exist. One is singu- 
larity subtraction. The form of the singular part of the solution is determined analyt- 
ically, and the solution is written as a sum of the singular solution and an unknown 
regular part. Only the regular part is treated numerically. The effect of the singular so- 
lution on the numerical one can be filtered down to the scales representable by the spa- 
tial resolution. The main advantage of this method is that it recovers the convergence 
of the scheme and at the same time approaches the exact solution. Recent applications 
of this method are to the driven cavity problem [34] and to injection of fluid into a 
cylindrical channel [35]. The results obtained are generally of high precision and often 
provide a benchmark for a particular problem. The main drawback is that it requires 
knowledge of the solution near the singular point. For the 2D driven cavity problem, 
the nature of the singularity was given by Dean and Montagnon [36] and Moffatt [37] 
for a Stokes flow. For most inertial (Navier-Stokes) flows, and for 3D flows Hills and 
Moffatt [38], the analytic form of the singular solution is unknown. Note that even 
when the velocity boundary conditions are continuous, lower-order singularities of 
purely geometric origin are present at the corners. We will return to this in section 3.2. 

Another approach is to derive a physically justified model which is no longer singu- 
lar. Methods from molecular dynamics reflect the microscopic nature of the fluid at 
the smallest scales but are very hard to adapt to problems containing both large and 
small scales. Several continuous (macroscopic) approaches have been proposed as a 
compromise between a continuous and a molecular description. These all introduce a 
spatially limited physical effect which effectively removes the singularity. In this cat- 
egory are methods based on variable slippage, as well as the surface viscosity or dy- 
namic surface tension applicable to free-surface problems. A comprehensive review of 
physically justified models, as well as other regularization techniques, is provided by 
Nguyen and co-workers [39, 40]. 

2.3 Time integration 

We use an implicit scheme for the linear diffusive terms while treating the nonlinear 
terms explicitly. Spectral methods require that the coefficients representing the solu- 
tion decay with their index or wavenumber. The nonlinear term in the Navier-Stokes 
equation can be seen as a generator and amplifier of high wavenumbers while the vis- 
cous term damps these high wavenumbers. The intensity of this damping depends 
on the particular time-integration scheme and on the way the Laplacian is evaluated 
and must be strong enough to oppose the effect of the nonlinear term. In our case, 
high-wavenumber modes are needed to represent both the thin boundary layers cre- 
ated near the rotating disks and the steep regularized boundary profile. Fortunately, 
these effects are most pronounced in the proximity of the boimdaries, where the ax- 
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ial Chebyshev and radial polynomial grid is finest. However, in the counter-rotating 
case, the central shear layer can also require high wavenumber modes in order to be 
well-represented . 



We use the jfirst-order backward Euler scheme for linear terms because it attenuates 
high frequencies faster than all other methods. Tests performed with the Crank-Nicolson 
method confirmed that for this scheme, nonlinear simulation was unstable even for 
quite small Reynolds numbers Re ~ 300. This behavior was also observed for the von 
Karman flow by Speetjens [8] and Lopez et al. [14] and by Marcus [41] in Taylor- 
Couette flow. The nonlinear term is treated by the second-order explicit Adams-Bashforth 
scheme, so that the backward Euler/ Adams-Bashforth time-integration scheme for the 
potentials and (p takes the following form: 

(7 - StRe-'A)Ahr^' = \r + y {^S^^ " SJ"^) ^ rhs^^ (2.12a) 
(7 - 3tRe-^A)AAh(l>''^^ = AA/j(/)" + y (^3SJ - Sj"^) = r/zsj (2.12b) 

where and SpSi are defined in (1.4). Equations (2.12) can be written as the nested 
system of equations: 

(l - 5tRe-^ a) = rhs^ (2.13a) 

Ah^ = (2.13b) 

- StRe-'^A^ g = rhs^ (2.13c) 

Aff^g (2.13d) 

Ah(p = U (2-13e) 



As explained in section 1, the boundary conditions imposed on (2.13) are Dirichlet 
conditions with boundary values calculated via the influence matrix in such a way as 
to satisfy the more complicated coupled boundary conditions given in (1.5). 

The maximal time step St depends on the Reynolds number. Typically starting from 
state u = requires a St which is 4-10 times smaller than that which can be used for 
evolving a fully developed state at the same Reynolds number. This is because the state 
u = is incompatible with the boundary condition (2.9a). In the first few iterations a 
boundary layer is created near the rotating cylinder lids, requiring higher spatial res- 
olution. This can be avoided by performing about 100 initial steps of the linear Stokes 
solver, i.e. with S,^ = = 0. 

In Table 1 we present the values of St and spatial resolutions typically used for per- 
forming nonlinear simulations for different values of Re. 
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Re 


configuration 


ot 


resolution (M x K x N) 


< 0(500) 


2D 


0.05 - 0.1 


1 X 32 X 16 


500 - 1000 


2D 


0.02 - 0.05 


1 X 64 X 32 


1000 - 3000 


2D 


0.01 - 0.02 


1 X 96 X 48 


3000 - 5000 


2D 


0.005 - 0.01 


1 X 128 X 64 


5000 - 10000 


2D 


0.001 - 0.0025 


1 X 180 X 90 


< 0(500) 


3D 


0.04 - 0.1 


8 X 64 X 32 


500 - 1000 


3D 


0.01 - 0.04 


16 X 80 X 40 


1000 - 3000 


3D 


0.025 - 0.01 


32 X 100 X 60 


3000 - 5000 


3D 


0.001 - 0.0025 


(64-96) X 128 X 80 



Table 1 

Typical values of timestep 5t and azimuthal (M), axial (X) and radial (N) resolution, for differ- 
ent configurations. Typical steepness of the regularization profile isS ^ 0.01 

2.4 Viscous terms 

We now describe the way in which the Helmholtz and Poisson problems with Dirichlet 
boundary conditions in (1.6) or (2.13) are solved. The azimuthal Fourier representation 
in (2.1) decomposes each 3D elliptic problem into a set of 2D problems, each of which 
is associated to a single azimuthal Fourier mode m. The reflection symmetry in z leads 
to further decoupling between the set of modes that are symmetric or antisymmetric in 
z. Each of the 2M resulting elliptic problems corresponds to a single azimuthal Fourier 
mode m and axial parity p G {s,a}, within which each equation corresponds to a value 
of {k,n), the indices of the axial and radial basis functions. The number of axial modes 
of each parity isK/2, and the number the number of radial basis functions correspond- 
ing to Fourier mode m is Nm = N — [y]; in the remainder of this section, we will take 
m = 0, so that Nm — N. 

In the full cylinder, regularity at the axis serves as one of the boundary conditions and 
is imposed by the use of the Qn radial basis in (2.1), leaving only the boundary condi- 
tion at r = 1 to be imposed. In the axial direction, the boundary conditions at the two 
disks can be recombined to yield one condition for each parity. Thus, one radial and 
one axial boundary condition remains to be imposed on each 2D problem. These are 
imposed via the t method [23, 24], so that the equations corresponding to the highest- 
wavenumber modes in each direction are replaced by the boundary conditions. 
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Gaussian elimination is especially economical for systems resulting from the spectral 
discretization of differential equations, whose structure is barely altered when bound- 
ary conditions are imposed via the t method, since recursion relations reduce the solu- 
tion time from quadratic to linear in the number of modes if the resulting systems are 
diagonally dominant. This can only be done in one direction however. In geometries 
with more than one non-periodic direction, the remaining directions must be treated 
by diagonalization. Here, we treat the axial direction by incorporating the boundary 
condition via Schur decomposition into the matrices representing 9^ for each parity 
and diagonalizing [24, 42], leading to decoupled problems for each axial eigenvalue 
Az- The operation count at each timestep, dominated by multiplication by the eigen- 
vector matrix, is quadratic in K/2 and linear in N. 

Thus, the 2D and 3D problems in (2.13) are all decomposed into a set of one-dimensional 
radial problems: 



which we will write in practice as 



hrdr -^ + ^ f-g (2.14) 



r^Hf= rdrrdr - + f = r^g (2-15) 

where m is the Fourier mode. The scalar A is or Az for the Poisson problems (2.13b), 
(2.13d) or (2.13e), or Az + Re/3t for the Helmholtz problems (2.13a) and (2.13c) (with 
the multiplicative factor —Re /St incorporated into g). 

With / and g represented in the polynomial radial basis (2.8), a recursion relation exists 
for r^H, as stated in section 2.1, i.e. r^H = R~^L with J^, L banded matrices. Thus each 
Helmholtz problem (2.15) can be replaced by 

Lf = Rr^g = Qg (2.16) 

For nonzero A, L is pentadiagonal and R is tridiagonal. Two obstacles must be sur- 
mounted before (2.16) can be solved in a time which is linear in the number of radial 
modes. First, L is not diagonally dominant, so that stable Gaussian elimination would 
require pivoting, destroying the banded structure. Second the radial boundary condi- 
tion must be included. A method which overcomes these difficulties was presented 
by Matsushima and Marcus [22]. Here, we will recast this method in terms of the 
Sherman-Morrison-Woodbury formula, which can be shown [21] to underly a large 
class of transformations between coupled and uncoupled systems. 

System (2.16) must be altered in two ways. First, since the largest element of L is located 
on the first super-digonal, permuting its rows leads to a matrix PL which is diagonally 



2, 
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dominant, but is no longer banded. Second, the boundary conditions must be imposed. 
The tau method replaces (2.15) by 



(2.17a) 
(2.17b) 



which effectively discards gj^ by adjusting it with the extra unknown t, introduced 
along with the boundary condition (2.17b). is the row vector which represents the 
discretized boundary condition and ej^ is the unit vector which selects the N^^ compo- 
nent. Multiplying by r^, using (2.16) and permuting rows leads to: 



r^Hf = R-^Lf = r^g + t^cnT 

Lf = Rr^g + Rr^e^r 
PLf^PQg + PQeNT 



(2.18) 



This is rewritten in matrix form as 



V 





—Qn,n\ 


/ \ 


/ iQ8)N \ 


^lo,i 


-Qlo,N 


/ 


-- iQg)io 




/ 


V ^ / 


\ ^ I 



(2.19) 



where the subscript lo refers to all indices lower than N. System (2.19) is solved by 
using the Sherman-Morrison-Woodbury formula 



(2.20) 



which relates the inverses of two matrices differing by a low-rank transformation vw^, 
in particular differing by a few rows or columns. We define A, v and for (2.19) as 
follows. The matrix in (2.19) is replaced by another matrix A, which is more easily 
inverted since it is block upper triangular: 



A 



( \ 

f 



( 





-Qn,n 


Llo,i 


-Qlo,N 





1 



/ 



/ iQg)N \ 




(2.21) 



where a is an arbitrary value whose order of magnitude is that of the dominant values 
of L and e\ is a unit row vector corresponding to the lowest radial wavenumber present 
for this m. The upper-left matrix in (2.21) is banded and diagonally dominant and so 
can be stably inverted without any pivoting. The matrices in (2.21) and (2.19) differ 
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only in their first and last rows, so their difference vw'^ is of rank two: 



/ Ljv,i - ael 




0\ / 1 \ 








vw^ (2.22) 



We define: 



Voiy 



I _ 







V 



, q = PQe^ 



I 



(2.23) 



Q/o,N 
\ / 

In addition to the ability to solve (2.21), the Sherman-Morrison-Woodbury formula 
(2.20) requires only the inversion of the following 2x2 matrix: 



/io\ 



-ae\ + e^L 







A" 




vol/ 

elL{L')-Hi {-ael + elL){L')-^qe^' 
B'^iL')-^ei BTiV)-^qeN 



(2.24) 



where we have used the fact that aeJ{L') ^ei — 1/a. 



2.5 Nonlinear terms 



To compute the nonlinear terms 

= Cz • V X S (2.25a) 
=-ez • V X V X S (2.25b) 

where 

1 

S = (u- V)u = -V(u-u)-uxa; (2.26) 

we use the pseudo-spectral method [23], in which fields are transformed into phys- 
ical space, the nonlinear terms are carried out via pointwise multiplication, and the 
results transformed back into spectral space. Computing the nonlinear term S in the 
rotational form — u x a; requires only 9 spectral <->physical transforms as compared to 
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the 15 transforms required by the convective form (u • V)u. The difference between 
them is annihilated by the curls taken in (2.25), so we will write S = — u x a;. 

The calculation of the nonlinear terms presents two difficulties. The first involves radial 
parity and appears when creating S. We have sought to use only scalar fields which can 
be represented by expansions of type (2.1) and constructed using radial operators such 
as (2.6) which preserve radial parity and can be implemented via recursion relations. 
The components of velocity and vorticity, defined in cylindrical coordinates using the 
toroidal and poroidal potentials as 



do not have this property. We therefore construct modified fields: 

u* = rurer + rueee + UzCz = {de^ + rdrz(p)^r + (90z^ - rdr^)e0 + {^h<P)^z (2.28a) 
CO* = rcorer + rcoeh + <^z^z = i^eUz - ^zUe)er + (dzU* - rdrUz)ee + {-^h^)^z 



whose components have the same parity as xp and ^ and so can be created and acted 
upon using the differential operators in (2.6). The modified fields u* and to* are trans- 
formed into physical space, where their cross product is taken to form 



The second difficulty appears when differentiating S in (2.25) and involves regularity. 
A vector function which is analytic at the origin must obey conditions analogous to 
(2.2), namely 



f,{r) = r\^-\r{r^) fe{r) = r^^Seir^) /z(r) = r^pzir^) (2.30) 



where pr, Pe and pz are polynomials. We require not only regularity at the origin of S, 
but also regularity of its curl and double curl. When S = — u x a; is calculated analyt- 
ically, this is in fact the case. However, the numerical transforms to and from physi- 
cal space introduce aliasing errors which destroy this property. Full dealiasing would 
multiply the time necessary for evaluating of the nonlinear term by a factor ~ 4.5. Mat- 
sushima and Marcus [22] suggest instead that all terms that could potentially suffer in 
spectral space from singular operations (like dividing by r) be evaluated in physical 
space (at collocation points excluding the coordinate origin) and transformed back to 
the spectral space using the radial transform, ensuring the correct polynomial order 
for a given Fourier mode. We have generalized this approach to the evaluation of Sip 




(2.28b) 



S* = rSrCr + rSeee + SzCz = S*er + SqCq + S^Cz = -u* x (v* 



(2.29) 
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and Sep. For each Fourier mode m, we write 

1 1 itn 

S^^ -{drvSe -imSr) -^{rdr - m)Sl - -^{S* + iS^) (2.31a) 

1 \ m 

(2.31b) 

The purpose of the decomposition in (2.31) can be explained as follows. Consider the 
first terms on the right-most-sides of (2.31), those which do not contain (S* + zS^). 
Because they do not generate terms of lower polynomial order, these terms preserve 
regularity and can be carried out in spectral space. In contrast, although (S* -|- fS^) 
should be divisible by r^, aliasing error in the multiplication (2.26) can generate terms 
of lower polynomial order. (More specifically, S* and iS^ are not individually divisible 
by and aliasing perturbs the cancellation in the sum S* + iSg.) This division by is 
therefore carried out in physical space. This increases the number of spectral ^-H'physical 
transformations only from 9 to 10. Further details can be found in [43]. 

2 . 6 Parallelization 

The separability of almost the entire algorithm (except for the nonlinear term) be- 
tween the Fourier modes makes parallelization of the code quite straightforward. Our 
code was parallelized using the MPI protocol which made it possible to run even very 
time-consuming three-dimensional simulations with resolutions such as M x KxN — 
128 X 160 X 90. Spectral methods are often considered to be poorly suited for paral- 
lelization as they require the exchange of all the data at each timestep of the simulation. 
In our code, all the necessary data exchange is done within two calls to the MPI_Alltoall 
MPI subroutine treating, in total, 10 three-dimensional fields at each time step. Even 
though this may seem to be a large operation, on the IBM Power4 architecture with 64 
processors we found that the time overhead per timestep due to the data exchange is 
counterbalanced by more efficient usage of the processor cache memory: each proces- 
sor of the parallelized code treats smaller data portions which can more easily fit into 
processor's fast internal memory (cache). We observed that the total CPU time used by 
the parallel code is often smaller than that used by the serial code treating the same 
problem, as shown in table 2. The efficiency of the parallel code depends, however, 
on the speed and the latency of the inter-processor network: the IBM Power4 architec- 
ture has particularly fast connection between the nodes which use the mixed model 
fast network/shared memory communication between processors. We conclude that for 
modem massively parallel computers, parallelization of a pseudo-spectral code does 
not necessarily degrade its efficiency but can actually enhance it. Additional technical 
information about the MPI parallelization of the code is given in [43]. 
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number of 


resolution 


CPU (sec)/ 


CPU (sec)/ 




processors 


MxKxN 


timestep 


(timestep x M) 


2D 


single processor 


1 X 96 X 48 


0.035 


0.035 


3D 


single processor 


32 X 96 X 48 


1.7 sec 


0.053 


3D 


32 processors 


32 X 96 X 48 


1.5 sec 


0.047 



Table 2 

CPU timings on IBM Power4, as the azimuthal resolution M and number of processors is var- 
ied. 

3 Tests and validation 

We now describe the ways in which we have validated the hydrodynamic code de- 
scribed here and in our companion paper [20]. We have obtained an exact polynomial 
solution to the nested Helmholtz-Poisson solver, which is the by far most complicated 
portion of the code; we present its form in the hopes it may prove useful to other re- 
searchers. For non-polynomial solutions, we have analyzed the effect of the corner sin- 
gularity on the error and verified the spectral convergence. We have tested the full non- 
linear time-dependent program by simulating 2D and 3D stationary and time-periodic 
rotor-stator flows which are well documented in the literature. 

3.1 Polynomial solutions 

There are no polynomial solutions which exactly satisfy the linear differential equa- 
tions (1.4) and boundary conditions (1.5) with which to compare a numerical solu- 
tion. We can, however, formulate a polynomial solution to a related problem: the time- 
discretized equations 

(I -StRe''^A)Ah^^rhsrp (3.1a) 
(J - StRe~'^A)AAh(l> = rhs^ (3.1b) 

relating the four fields tp, (p, rhsip, rhs<p and subject to boundary conditions (1.5). Un- 
der these conditions, no error is introduced in imposing the boundary conditions via 
the T method (see section 2.4). Comparison between this analytic solution and a nu- 
merical solution provides a test of the linear Helmholtz/Poisson solver, including the 
enforcement of the boundary conditions via the influence matrix method, which is in- 
dependent of the evaluation of the nonlinear terms, the temporal integration and the 
spatial resolution. 
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Using a symbolic algorithm implemented in Maple, we have calculated polynomial 
solutions which contain several radial and axial basis functions (wavenumbers) and 
which correspond to a realistic projfile cv±, such as o;^"'^ = ±(1 — r^). Polynomial so- 
lutions were calculated for Fourier modes ranging from m — to m — 5. For m — 0, we 
sought a solution containing two recirculation rolls separated by the mid-plane z — 0, 
leading to the potentials, right-hand-sides, and velocities: 

^pP°^y{r,z) = ^z(-30z2 + 33z4 + 5)r^ - -^z(z - l)(z + l){5z^ - l)r^ - \z^r^ 
64 48 2 

(3.2a) 

r}is^°^^(j,z) = -z(-690z2 + SSz* + 185)r^ + -z{-1970z^ + 425 + 1609z4)r4 (3.2b) 
- 60z(z - 1) (z + 1) (5z2 - - 40z^ + 2z^ (3.2c) 
(pPoly{r,z) = -\{r- ifir + - + (3.2d) 

rhs^^"^^{r,z) = -72z{5z^ - 33)r* - 96z(3z^ + 108 - ISlz^y + 1248z5 (3.2e) 
+ 4344Z - 6456z^ (3.2f) 
u^/y = -3r(z - l)(z + l)(5z2 - l)(r - i fir + i f (3.2g) 

-I 

ul°^\r,z) = --zr{r - l)(r + l)(-30rV + 5r^ + 33rV + 8rV + 8z^) (3.2h) 
uf^^ 6z(z-l)2(z + l)2(r-l)(r + l)(3r2-l) (3.2i) 

For all polynomial solutions tested, we obtained numerical solutions for which the 
maximum relative errors in tp, (p and in satisfaction of the boundary conditions, as 
well as the equations, are typically O (10-^^), i.e. machine precision, and never exceed 
0(10-12). 



3.2 Non-polynomial solutions: error analysis 



If the right-hand-sides rhs^ and rhstp are not polynomials that can be exactly repre- 
sented within the spatial discretization, then the equations cannot be exactly satisfied. 
Note that the tau error introduced in favor of satisfying the boundary conditions at 
each level of the nested system of equations is necessarily propagated to the next level. 
That is, tp"~^^ is not an exact solution to the Poisson problem (2.13b) and in addition, 
the right-hand-side is not an exact solution to the Helmholtz problem (2.13a). This 
implies that the error is not isolated in the equations corresponding to the highest 
wavenumbers, but is distributed among all the equations. However, we will see be- 
low that non-satisfaction of the equations for low wavenumbers does not have severe 
consequences. 
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In figures 6 and 7, we display the relative error 

\EAjitp — rhstj,\ \EAAi^(p — rhS(p\ 

sup I r^s^ I ' sup \rhS(p\ 

in physical space {r,z) for the m = and m = 2 modes. In this subsection, as in the 
previous one, we restrict ourselves to the Stokes problem, whose errors behave like 
those of the Navier-Stokes equations. As is usually the case, the error is concentrated 
in the neighborhood of the boundaries, decaying rapidly away from them. 




Fig. 6. Error of satisfaction of equations for ni = after 100 timesteps of nested Stokes solver. 
Upper surfaces show the error for the Stokes problem, lower surfaces for the polynomial solu- 
tion (3.2). Resolution used: K = 64,N = 32. Left: e"'=°. Right: e^=^. 




Fig. 7. Error of satisfaction of equations for m = 2, after 100 timesteps of nested Stokes solver. 
Upper surfaces show the error for the Stokes problem, lower surfaces for a polynomial solution. 
Resolution used: K = 6A, N = 32. Left: ef=^. Right: ef^=2. 

These errors represent quite severe criteria, since these measure the satisfaction of the 
curl and double curl of the original equations; the error in satisfying the Stokes equa- 
tions themselves is considerably lower, on the order of e^p divided by 0{K,N) and 
divided by 0{K^,N^) near the boundary. This estimate relies on the fact that the spec- 
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trum of the error in satisfying the Stokes equations is approximately uniform, so that 
its derivatives are dominated by its high wavenumber terms. 

The error depends significantly on the right-hand-sides rhsjp and r/zs^. For arbitrary 
right-hand-sides, the boundary conditions can be very constraining and can lead to 
a nearly singular solution suffering from spurious oscillations, with an error near the 
boundary that is 0(1). However, are considerably smaller when the right-hand- 

sides are calculated from the solution at the previous timestep, especially for the Stokes 
problem for which the nonlinear term is zero. 

For the axisymmetric modes shown in figure 6, e^p is O(10~^) on the boundaries and 
0(10^^*^) for the internal points. This is much less than e^, which reaches 0(1) at 
the cylinder comers. To understand this, we recall that in the axisymmetric case, the 
toroidal flow described by ip is azimuthal and the poloidal flow described by cp is in the 
(r,z) plane. The azimuthal flow described by xp follows smooth paths, while the flow 
described by (p must abruptly change direction near the comers. This poloidal flow in 
fact resembles the analytic asymptotic solution derived by Moffatt [37] for 2D flow in 
a rectangular container, which is weakly singular in that its vorticity behaves like p^-^^ 
(for a small distance p from the comer). This in turn implies that the axial component 
of the Laplacian of the Stokes equation measured by diverges for exact solutions to 
the continuous 2D Stokes problem, in contrast with the numerically computed solution 
which is forced to be finite and regular. Satisfaction of the Stokes equation itself follows 
from the satisfaction of the boundary condition which our code imposes to precision 
0(10-14). 

For the non-axisymmetric modes shown in figure 7, typical errors at the corners are 
0(0.01) for ip and O(O.l) for (p and the same corner singularity is observed for both. 
This is because the 3D Stokes solutions are also weakly singular at the corners Hills and 
Moffatt [38] and ip and (p are coupled for m 7^ 0. The relative error decreases rapidly 
away from the boundaries: it is 0(10^^) only two gridpoints away and reaches reaches 
O(IO^) for the interior points. 

The accuracy could be further enhanced - or confined to certain modes - by completing 
the method with the T-correction [21], which takes into account the high-wavenumber 
residuals resulting from the imposition of the boundary conditions for each Helmholtz 
or Poisson problem. 

3.3 Spectral convergence 

The important indicator of spatial convergence for spectral methods is the decay rate of 
high-wavenumber coefficients in the solution fields. For a well-behaved solver, in the 
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Fig. 8. Spectral coefficients for m = 0, after 100 timesteps of the nonlir\ear Navier-Stokes solver 
{Re = 750, St = 0.01). Resolution used: K = 64, N = 32. Left: ^5/"'=° Right: pftp'^^. 

absence of volume and boundary singularities, the magnitude of spectral coefficients 
should decay rapidly with wavenumber. For a laminar flow this decay rate should be 
exponential, but the presence of thin boundary layers can significantly influence the 
convergence. Figure 8 shows the spectral convergence for a full Navier-Stokes simula- 
tion at Re = 750 at T = lOOSt = 1 from initial conditions of uq = (x)j^{r)2z/h,Uz = Ur = 0. 
We show the {r,z) spectral coefficients for Fourier mode m = 0; those for m 7^ are 
similar. The convergence of the spectra can be qualified as quasi-exponential, mean- 
ing that the high-wavenumber spectral coefficients seem not to decrease below a level 
can almost certainly be attributed to the singular character of the solu- 
tion to the Stokes equation near the cylinder comers. In support of this explanation, 
we note that the spectrum obtained after 100 timesteps of the linear (Stokes) solver 
behaves similarly, except for xp"^^^, which displays true explonential decay, down to 
levels of 10~^^ for the same resolution. 

3.4 Axisymmetric rotor-stator configuration 

The code was first tested on the well-documented axisymmetric rotor-stator configura- 
tion with aspect ratio h = 2. The first test is the reproduction of the characteristic steady 
state for Re ~ 1850 where the flow exhibits two recirculation bubbles (one large and the 
other much smaller) situated approximately at (r = 0, z = 1/2) and (r = 0, z = 0). The 
contour plot of the Stokes streamfunctions defined as 

(7(r,z) = —rdr'^(r,z) ^{r,z) = —rdr(p{r,z) (3.3) 

presented on figure 10 matches that presented by Daube [6] (figure 9) and is similar 
to those of Lopez and Shen [7] obtained for the slightly larger aspect ratio h = 2.5. 
Quantitative agreement between our results and the previous calculations by Daube 
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Fig. 9. Contours of Stokes Fig. 10. Rotor-stator configuration for Re — 1850 and h = 2. 
poloidal streamfunction t] Left: poloidal streamf unction rj. Right: toroidal streamfunc- 
from [6] for Re = 1850, h = 2. tion cr. 



[6] and Lugt and Abboud [3] is established by comparing the profiles of axial velocity 
Uz on the cylinder axis (see figure 12). This test shows excellent agreement between 
our results obtained using the poloidal-toroidal formulation (fig 12b) and the velocity- 
vorticity formulation (fig. 12a). 

It was observed experimentally by Escudier [44] and numerically by Daube and Sorensen 
[4], Lopez [5], Daube [6], Gelfgat et al. [10] and Speetjens [8] that this flow undergoes a 
Hopf bifurcation toward a flow oscillating with a frequency approximately 0.25 times 
the rotation frequency. This transition occurs at Re near 2600 with a period of T = 26.55. 
These values are in the ranges previously found for this configuration (see table 3); de- 
viations can probably be attributed to the differences in the treatment of the boundary 
conditions (regularization). The simulation was performed with St = 0.01 and with 
high spatial resolution K x N = 140 x 70 in order to well represent the sharp regu- 
larization profile corresponding to 3 = 0.06 imposed on the lateral boundary r = 1 (see 
section 2.2) as proposed by Lopez and Shen [7] and also used by Speetjens [8]. The time 
evolution ug(r — 0.5, z — 0,t), along with the normalized power spectrum, are shown 
on figures 13 and 14. 
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Method Reference 



u — 


CO 


Daube [6] 


25.52 


V- 


CO 


Daube [6] 


25.84 


u — 


CO 


Speetjens [8] 


26.61 


u — 


V 


Gelfgat et al. [10] 


f^26.7 


^- 


<p 


this work 


26.55 



Table 3: Oscillation period in rotor-stator 
configuration with h = 2 and Re = 2800. 




Fig. 11: Rotor-stator configuration with 
Rec ~ 2150 and h — 3.5. Isosurface Mz ~ 
of axial velocity field after subtraction of 
axisymmetric component. 
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Fig. 12. Profile of at r = for rotor-stator configuration at Re = 1850, h = 2. a) Profile from 
Daube [6]: results obtained using rj — co {*) and u — a; (o) are superposed, b) Profile att = 3000 
obtained from the present code, using poloidal-toroidal decomposition ^ — (p- 



3.5 First 3D instability 



We have tested the non-axisymmetric aspects of the code using another rotor-stator 
configuration. For the aspect ratio oih = 3.5 we found that the first bifurcating mode 
has wavenumber m = 3 and critical Reynolds number Rcc = 2116. This result is in good 
agreement with Gelfgat et al. [10], where the critical Reynolds number was estimated 
at Re — 2131. A characteristic spiral analogous to that visualized by Lopez et al. [14] 
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t 

Fig. 13. Saturated state of time evolution of 
UQ{r = 0.5, z = 0) for the rotor-stator configu- 
ration at Re = 2800, h = 2. 




Fig. 14. Density power spectrum of 
ug{r = 0.5, z = 0) from figure 13. 



0.1 
0.08 
0.06 
0.04 
0.02 
0, 




1000 2000 3000 4000 5000 
t 




1000 2000 3000 4000 5000 
t 



Fig. 15. Time history of ug{r = 0.5, z = 0) for the rotor-stator configuration at Re = 2800, h = 2. 
Left: from Speetjens [8]. Right: current work. 



for the same configuration is represented on figure 3.4. 



4 Conclusion 



Motivated by a need for a numerical tool adequate for investigating cylindrical von 
Karman flow, we have written a spectral code which solves the Navier-Stokes equa- 
tion for an incompressible fluid in a finite cylinder. This task encompasses a number 
of algorithmic challenges, which we list in order of decreasing difficulty. The first is 
to impose incompressibility, a goal which we achieved by using the poloidal-toroidal 
decomposition. In a finite cylinder, this formulation leads to differential equations of 
higher order and coupled boundary conditions [18, 19]. In a companion paper [20], we 
described the way in which we used the influence matrix technique to transform these 
equations and boundary conditions into an equivalent set of decoupled Helmholtz or 
Poisson problems, each with Dirichlet boundary conditions. 
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The second challenge, which has been the main focus of this article, is the treatment of 
the cylindrical axis. The singularity engendered by the use of cylindrical coordinates 
is only apparent and should not be transmitted to the fields. We have dealt with the 
axis singularity by using a basis of radial polynomials developed by Matsushima and 
Marcus [22] which are analytic at the axis and which have properties similar to those 
of the Legendre polynomials. We have extended [22] in several ways, (i) The radial 
basis was developed for two dimensions, i.e.polar coordinates. Here, we have used it 
to represent fields in a three-dimensional cylinder of finite (non-periodic) axial length. 
This fairly straightforward extension was implemented by diagonalizing the differen- 
tial operators in the axial direction, (ii) Numerical inaccuracy in the pseudo-spectral 
evaluation of the nonlinear term can generate terms which are not analytic. We have 
generalized to the high-order equations resulting from the poloidal-toroidal decompo- 
sition the procedure developed by Matsushima and Marcus [22] to avoid this problem. 
Care must be taken to preserve order and parity at each stage of the calculation, (iii) 
Differential operators expressed in this basis can, as in most such cases [31], be refor- 
mulated as recursion relations, which can be used to reduce the time for action or in- 
version of differential operators. Using the Sherman-Morrison- Woodbury formula, we 
have formalized the stable algorithm given in [22] for solving Poisson and Helmholtz 
problems in a time proportional to the number or radial gridpoints or modes. 

The third challenge is the genuine singularity at the corners of the domain, where the 
disks and the cylinder which bound the domain meet. Finite difference and finite el- 
ement codes intrinsically smooth the singularity; in contrast, spectral expansions at- 
tempt to converge to the discontinuity, leading to spurious oscillations. In our code, 
we replaced the discontinuous boundary conditions by a profile on the disks which is 
steep but continuous as the corner is approached. The geometrical singularity remains, 
but it is weak and does not prevent spectral convergence. 

Tests performed for analytic polynomial solutions to the Helmholtz problem with an 
appropriate right-hand side showed that the solver reproduces exact solutions to nearly 
machine precision. For a non-polynomial solution, the solver displays exponential con- 
vergence of spectral coefficients of the solution. The potential equations (correspond- 
ing to the curl and double curl of the Navier-Stokes equations) are satisfied, with an 
error of only O(10~^°) for the interior points. The numerical code was parallelized 
using the MPI protocol. This made it possible to simulate nearly turbulent flow for 
Re = 5000 with a spatial resolution of 128 x 160 x 90. 

Finally, we have validated the hydrodynamic code by testing it against well-documented 
problems in the literature, demonstrating the feasibility of calculating solutions to the 
time-dependent Navier-Stokes equations in a finite cylindrical geometry which are 
both analytic and divergence-free. 
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